N=1;        % Truncation order
M=1;        % Number of polynomial terms (M=3 is sufficient).

index_1= [  1 ] ;

for ind1= 1:length(index_1)
    
    index_2  = index_1(1,ind1);
    
period_1=linspace(  1 ,  1 ,  1 );     %% 角度


for ww1=1:length(period_1)
      pp=period_1(1,ww1);   % wavelength (um)


angle_1=linspace( 0,  0  ,  1  );     %% 角度



for ss1=1:length(angle_1)

angle_spr=   0  ;

      
 Audata1=load('TiTi2.txt');          %% 金属材料折射率   

lamb_1=linspace(   0.4  ,  2   , 300 );     %% 入射波长 

for llm1=1:length(lamb_1)         %% 可以并行运算
      lam01=lamb_1(1,llm1);   % wavelength (um) 
  a=lam01*1000;



n_index1 = Audata1(llm1,2)+Audata1(llm1,3)*1i;   %% 计算的标志金属折射率的虚部要负值 所以修改这样
      
        
p=pp;   %% 周期
lambda = lam01/p  ;     % Free space wavelength normalized to grating period
n1 = 1      ;         % Refractive index of incident medium
n3 = 1  ;     % Refractive index of transmission medium
eps1=n1*n1;     %Permittivity of incident medium
eps3=n3*n3;    %Permittivity of transmission medium
alpha=angle_spr*pi/180+1e-6;% Angle of incidence (in radian)

 t6=     2 ;





%tttttt= t1  ;



d_1=  t1/p  ;
d_2=   t2/p  ;
d_3=  t3/p  ;
d_4=   t4/p  ;
d_5 =  t5/p ;
d_6=  t6/p ;

fill=0.5;
w=   p*fill   ;
xx1=  w/p ; 


n_SiO2=  1.46*1.46 ;
n_Al2O3=  1.73*1.73 ;
n_TiO2=  2.41 *2.41  ;



n_Ti = (n_index1)*(n_index1)  ; 


xt=  [ xx1 xx1 xx1 xx1 xx1 xx1]  ;
epst= [   n_SiO2 n_SiO2
    n_Al2O3 n_Al2O3
    n_TiO2 n_TiO2
    n_Ti n_Ti
    n_SiO2 n_SiO2
    n_Ti n_Ti

  
          ].'   ;    
      
dl=[ d_1 d_2 d_3 d_4 d_5,d_6];

 
% -------------------
% initial computations
% -------------------
k0=2*pi/lambda;
k1=k0*(eps1^.5);
k3=k0*(eps3^.5);
ux=cos(alpha);
uz=-sin(alpha);
kx = k1 * sin(alpha) - (-N:N) * 2 * pi;
kz1 = sqrt(k1^2-kx.^2);
kz3 = sqrt(k3^2-kx.^2);
for j=1:2*N+1
    if real(kz1(j))>0
       kz1(j)=-kz1(j);
    end
    if imag(kz1(j))<0
       kz1(j)=-kz1(j);
    end
    if real(kz3(j))<0
       kz3(j)=-kz3(j);
    end
    if imag(kz3(j))>0
       kz3(j)=-kz3(j);
    end
end
Kz1=diag(kz1);
Kz3=diag(kz3);
% ------------------
 
Nlayer = length(dl);
for lcounter=1:Nlayer
    d=dl(lcounter);
    if max(epst(:,lcounter))==min(epst(:,lcounter))
        [r11 r12 r21 r22]=RmatrixTM_hom(epst(1,lcounter),d,lambda,alpha,eps1,N);
    else
        [r11 r12 r21 r22]=RmatrixTM_lam(xt(:,lcounter),epst(:,lcounter),d,lambda,alpha,eps1,N,M);
    end
 
%   -------------------------------
%   R matrix Propagation algorithm
%   -------------------------------
    if lcounter==1
         R11=r11;
         R12=r12;
         R21=r21;
         R22=r22;
    else
         Z=inv(r22-R11);
         Rx11=r11-r12*Z*r21;
         Rx12=r12*Z*R12;
         Rx21=-R21*Z*r21;
         Rx22=R22+R21*Z*R12;
         R11=Rx11;R12=Rx12;R21=Rx21;R22=Rx22;
    end
%   -------------------------------
end
 
% ---------------------------------------------------------------
% Finding reflection and transmission coefficients using R matrix
% ---------------------------------------------------------------
GG3=diag((kx.^2+kz3.^2)./kz3);
GG1=diag((kx.^2+kz1.^2)./kz1);
 
 
I = eye(2*N+1);
G11=I-R22*GG1/k0;
G12=-R21*GG3/k0;
G21=-R12*GG1/k0;
G22=I-R11*GG3/k0;
G=[G11 G12;G21 G22];
 
deltai0=zeros(2*N+1,1);
deltai0(N+1)=1;
b=0;
uprim=(k1*cos(alpha)*ux-k1*sin(alpha)*uz);
b(1:2*N+1)=(-ux*I+uprim*R22/k0)*deltai0;
b(2*N+2:4*N+2)=(uprim*R12/k0)*deltai0;
RT=G\b.';
 
Rx=RT(1:2*N+1).';
Tx=RT(2*N+2:4*N+2).';
 
Rz = - kx .* Rx ./ kz1;
Tz = - kx .* Tx ./ kz3;
% ---------------------------------------------------------------
    
% --------------------------------------
% Evaluation of Diffraction efficiencies
% --------------------------------------
DE1=zeros(1,2*N+1);
DE3=zeros(1,2*N+1);
for j=-N:N
    DE1(j+N+1)=-real(kz1(j+N+1)/k1/cos(alpha))*(abs(Rx(j+N+1))^2+abs(Rz(j+N+1))^2);
    DE3(j+N+1)=real(kz3(j+N+1)/k1/cos(alpha))*(abs(Tx(j+N+1))^2+abs(Tz(j+N+1))^2);
end


RR(1,llm1)=sum(DE1);      %%  计算反射率
AA(1,llm1)= 1- RR(1,llm1); 

fprintf('\n 菜园疯狂跑了 %d 公里',llm1 )
end
end
end
end